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Abstract: The present article proposes an approach to analyze the photometric properties of the surface materials 
from multi-angle observations acquired by the Compact Reconnaissance Imaging Spectrometer for Mars (CRISM) on- 
board the Mars Reconnaissance Orbiter. We estimate photometric parameters using Hapke's model in a Bayesian inversion 
framework. This work also represents a validation of the atmospheric correction provided by the Multi-angle Approach for 
Retrieval of Surface Reflectance from CRISM Observations (MARS-ReCO) proposed in the companion article. The latter 
algorithm retrieves photometric curves of surface materials in reflectance units after removing the aerosol contribution. 
This validation is done by comparing the estimated photometric parameters to those obtained from in situ measurements 
by Panoramic Camera instrument at the Mars Exploration Rover (MER)-Spirit and MER-Opportunity landing sites. 
Consistent photometric parameters with those from in situ measurements are found, demonstrating that MARS-ReCO 
gives access to accurate surface reflectance. Moreover the assumption of a non-Lambertian surface as included in MARS- 
ReCO is shown to be significantly more precise to estimate surface photometric properties from space in comparison to 
methods based on a Lambertian surface assumption. In the future, the presented method will allow us to map from orbit 
the surface bidirectional refiectance and the related photometric parameters in order to characterize the Martian surface. 

1 Introduction 

Refiectance of planetary surfaces is tightly controlled by the composition of the present materials but also by their 
granularity, the internal heterogeneities, porosity, and roughness. The refiectance can be characterized by measurements 
at different wavelengths, viewing geometries (emergence direction) and solar illuminations (incidence direction). Such 
investigations have been conducted for Mars using telescopes, instruments on-board spacecrafts and rovers. A summary 
of these studies is available in the Johnson et al. ( |2008 1 review chapter. One can find studies related to: the Viking 



1999[ ), the Hubble Space Telesco pe ^Bell III 



Landers {Guinness et al. 19971, the Pathfinder Lander (Johnson et al. 

et al. , 1999), the PANoramic CAMera (Pancam) instrument on-board Mars Exploration Rovers (MER) ( Johnson et 
al. |2006a^b) , the Observatoire pour la Mineralogie, I'Eau, les Glaces et lActivite (OMEGA) instrument on-board Mars 



Express (MEx) ( Pinet et al. 2005 1 and the High Resolution Stereo Camera (HRSC) instrument on-board MEx ^Jehl et al. 
2008). Recently, Shaw et al. (2012) derived maps of mm- to cm- scale surface roughness at MER-Opportunity landing site 



by using multi-angle hyperspectral imager spectrometer called Compact Reconnaissance Imaging Spectrometer for Mars 
(CRISM) on-board Mars Reconnaissance Orbiter (MRO). 

In order to derive compositional and structural information from refiectance measurements, physical models describing 
the interaction of light with natural media are needed. Chandrasekhar ( [i960 ) proposed the radiative transfer equation 
describing the loss and gains of multidirectional streams of radiative energy within media considered as continuously 
absorbing and scattering where grains are separated by a distance greater than the wavelength (e.g. atmospheres). In 
the case of a dense medium (e.g., surfaces), two different solutions are developed. First, solutions based on Monte Carlo 
ray tracing methods handled the medium complexity (e.g. 



Grynko and Shkuratov 



2007). Unfortunately, this direct 
approach requires large computing times and large parameter space, limiting the inversion. Second, solutions based on 
an empirical or semi-empirical approach were proposed by adapting the radiative transfer equation to granular media 
(e.g., 



Hapke 1981a.. Hapke 1981h Hapke 



1986 


1993 


2002 


Shkuratov and Starukhina 


1999 


Doute and Schmitt 


1998 



These techniques are more relevant for inversion. Several parameters characterize natural surfaces such roughness and 
compaction, while other parameters characterize an average grain, such the single scattering albedo or the phase function. 

Previous photometric studies suggest that variations in scattering properties are controlled by local processes. For 
example, photometric variations observed by Pancam at Columbia Hill and the cratered plains of Gusev Crater are mainly 
caused by aeolian and impact cratering processes ( Johnson et al. 2006a ) . These conclusions encourage us to expand the 



estimation of photometric properties from in situ observations to the entire planet using orbital data to go further in the 
interpretations. 

Observations acquired from space, however, require the correction for atmospheric contribution (i.e., gases and aerosols) 
in the remotely sensed signal prior to the estimation of the bidirectional reflectance of the surface materials. Previous orbital 
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Figure 1: Detailed scheme of the estimation of surface photometric properties. The blue blocks represent the initial FRT 
CRISM observations (from i=l to n). The green blocks represent the aerosol optical thickness retrieval and the correction 
for atmospheric contribution for the determination of the surface bidirectional reflectance, which is carried out by the 
methodology described in the companion paper ( Ceamanos et al. 2013 1. The red blocks correspond to the work presented 



in this article, that is the estimation of surface photometric properties. 



photometric studies of Martian surfaces were conducted without atmospheric correction but using the lowest aerosols 
content observations (e.g., Jehl et al. 2008 Pinet et al. 2005[). Usin g data acquired by the multi-angle hyperspectral 



imaging spectrometer called CRISM on-board MRO {Murchie et al. 20071, our objective is to estimate accurate (i) 



surface bidirectional reflectance of the surface of Mars and (ii) photometric parameters associated with the materials. 
Ceamanos et al. ( 2013[ ) presents a method referred to as Multi-angle Approach for Retrieval of Surface Reflectance from 
CRISM Observations (MARS-ReCO) . This original technique takes advantage of the multi-angular capabilities of CRISM 
to determine the bidirectional reflectance of the Martian surface. This is done through the atmospheric correction of the 
signal sensed at the top of atmosphere (TOA). We propose an approach to analyze the photometric parameters of the 
surface materials in terms of structural information by inverting Hapke's photometric model in a Bayesian framework, as 
discussed below. The validation of the methods proposed in this work and in the companion article ( Ceamanos et al. 



2013 ) is performed by comparing the estimated photometric parameters to those obtained from in situ measurements by 
Pancam instrument at the MER-Spirit and MER-Opportunity landing sites (respectively at Gusev Crater and Meridiani 
Planum) ( Johnson et al. 2006a|b ) . 

This article is organized as follows. First, the methodology to obtain photometric surface parameters is described in 
Section [2j Second, the estimated photometric parameters are presented in Section [3j Third, results are compared to 
experimental studies, independent orbital measurements and in situ measurements in Section [4] The signiflcance of the 
photometric results shall be discussed in Section [s] 



2 Methodology 

This article and its companion take advantage of the multi-angular capabilities of the CRISM instrument to correct 
for atmospheric contribution in order to estimate the surface bidirectional reflectance {Ceamanos et al. 20131 and to 
determine the surface photometric parameters (this work). The approach presented in this article includes the following 
steps: (i) the selection of appropriate CRISM observations at both MER landing sites for the photometric study, (ii) 
the determination of the surface bidirectional reflectance by correcting for aerosol contributions, (iii) the combination of 
several CRISM observations for a better sampling of the surface bidirectional reflectance, and (iv) the estimation of the 
associated surface photometric parameters. The detailed scheme of the procedure is illustrated in Figure [l] One should 
note that, in order to test the performance of the method presented throughout this article and its companion paper, 
the study is only conducted at one wavelength and for some spatial pixels. We choose to work at 750 nm where (i) the 
contribution of gases is minimal and thus the retrieval of photometric properties is likely to be more accurate and (ii) in 
situ photometric measurements from Pancam instrument are available for the comparison to the estimated photometric 
parameters. 
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Table 1: Selected CRISM observations focused on the Spirit and Opportunity landing sites, respectively at Gusev Crater 
and Meridiani Planum. Ls, stands for the solar longitude, 6q is the incidence angle, g is the phase angle range. AOTminerai 
stands for the mineral Aerosol Optical Thickness at 1 fim from Wolff's estimates ( Wolff et al. 2009 1 (personal communica- 
tion of Michael Wolff) and AOT„ater is the water Aerosol Optical Thickness at 320 nm by MARs Color Imager instrument 
(personal communication of Michael Wolff) . 



Gusov Crater (MER-Spirit) 



Meridiani Planum (MER-Opportunity) 





FR,T3192 


FRTSCEl 


FRTCDA5 


FRT95B8 


FRT334D 


FRTB6B5 


Acquisition date 


2006-11-22 


2007-12-17 


2008-10-07 


2008-01-11 


2006-11-30 


2008-07-08 


Ls (degree) 


139.138 


4.040 


138.333 


16.223 


142.975 


96 


00 (degree) 


60.4 


40.02 


62.8 


39.3 


55.4 


56.4 


g (degree) 


~56-112 


~41-90 


~46-106 


~41-86 


~41-106 


~40-106 


AOT„,i„orai (1/Jin) 


0.33±0.04 


0.98±0.15 


0.32±0.04 


0.56±0.09 


0.35±0.04 


0.35±0.04 


AOT„at„ (320 nm) 


0.08±0.03 


0.07±0.03 


0.03±0.03 


0.12±0.05 


0.12±0.03 


0.14±0.03 



2.1 CRISM data sets 

2.1.1 The CRISM instrument and targeted observations 

The CRISM instrument on-board MRO is a visible and infrared hyperspectral imager (i.e. 362 to 3920 nm at 6.55 
nm/channel) that operates from a sun-synchronous, near-circular (255 x 320 km altitude), near-polar orbit since November 
2006. The appropriate mode to estimate surface spectrophotometric properties is the so-called targeted mode providing 
Full Resolution Targeted (FRT) observations consisting of a sequence of eleven hyperspectral images from a single region 
acquired at different emission angles. The solar incidence angle is almost constant during the MRO flyby of a targeted 
observation. A typical targeted sequence is composed of a nadir image (~10xlO km) at high spatial resolution (15-19 
m/pixel) and ten off nadir images with a lOx-spatial binning (resulting in a resolution of 150-200 m/pixel) taken before 
and after the nadir image. The latter sequence constitutes the so-called Emission Phase Function (EPF) sequence. The 



pointing of CRISM can rotate (gimbal) ± 60° (Murchie et al.\ 20071. 



2.1.2 Selection of targeted observations 



As explained in the companion paper ( Ceamanos et al. 2013[ ), the accuracy of the surface reflectance provided by MARS- 



ReCO when dealing with a single targeted observation highly depends on the combination of a moderate atmospheric 
opacity (i.e. aerosol optical thickness less than or equal to 2), reasonable illumination conditions (i,e. incidence angle less 
than or equal 60°), an appropriate phase domain (i.e. significant difference between the available maximum and minimum 
phase angles up to 40°) and on the number and diversity of angular measurements. The combination of several targeted 
observations, since it enables a better sampling of the bidirectional reflectance, could therefore significantly improve the 



refiectance estimation as it provides more regular angular sampling of the surface target. Pinet et al. (2005) and Jehl 



et al. ( 2008 1 proved the beneflts of using different spaceborne observations under varied illumination conditions (OMEGA 
and HRSC). The principal requirement to combine targeted observations is the absence of seasonal changes among the 
selected observations. 

Several CRISM observations have been acquired over the MER landing sites since the beginning of the mission. In 
particular, up to sixteen and ten CRISM full targeted observations (FRT) are available in the MER-Spirit and MER- 
Opportunity landing sites, respectively. In this article, we select CRISM observations according to several criteria: (i) the 
quality of overlap among the observations (above 70%), (ii) the variation of the solar incidence angle implying a widening 
of the phase angle domain (note that only the variation of the seasonal solar longitude (Ls) can provide different incidence 
angles due to the sun-synchronous orbit of MRO, and (iii) the absence of surface changes (e.g., seasonal phenomena) as 
they can jeopardize the determination of the surface photometric properties. Taking into account these criteria, three 
CRISM observations acquired over Gusev Crater (i.e., FRT3192, FRTSCEl and FRTCDA5) and over Meridiani Planum 
(i.e., FRT95B8, FRT334D and FRTB6B5) are selected, respectively. We note that the selected observations have quite 
different phase angle ranges as shown in Table [T] 

Targeted observations are archived in the Planetary Data System (PDS) and are composed of: (i) Targeted Reduced 
Data Records (TRDR), which store the calibrated data in units of I/F (radiance factor - RADF, see Table|2|, the ratio of 
measured intensity to solar flux, and (ii) Derived Data Records (DDR), which store the ancillary data such as the spatial 
coordinates (latitude and longitude) and the geometric configurations of each pixel by means of the incidence, emission 
and phase angles. In the present study, CRISM products are being released with the TRDR2 version of CRISM calibration 
(TRR2 for brevity). 
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Table 2: Photometric units derived from the bidirectional reflectance r. (/: intensity, F: solar flux, : 9q: incidence angle, 
9: emergence angle, g: phase angle). The CRISM observations are released in RADF unit and the SPC cubes in BRF 



units. 



Unit Symbol Name Expression 

r bidirectional reflectance r{9Q^0,g) — 

BB.DF bidirectional reflectance distribution function BRDF — r{9o^6,g)/ cos (^^o) 

eflectance sr 

RADF Radiance factor RADF = tt X r(0o, 9,g) = j, 

BRF bidirectional reflectance factor p — BRF — tt X r{9o^ 0^ g)/ cos (Oq) 



2.1.3 SPC cubes: Integrated multi-angle product 



To facilitate the access to the multi-angular information pertaining to each terrain unit, the eleven hyperspectral images 
corresponding to a single targeted observation were spatially rearranged into data set named SPC (spectro-photometric 
curve) cube (see Ceamanos et al. ( 2013[ | for more detail). The SPC cube is composed of: (i) x dimension (horizontal) 
corresponding to angular configurations (up to eleven) grouping all refiectance values at different geometric views, (ii) y 
dimension (vertical) corresponding to the spatial coordinate defining a super-pixel. The super-pixels are rearranged as 
function of the number of available angular configurations and, (iii) z dimension corresponding to the spectral sampling. 
The photometric curves, stored in a SPC cube are in units of RADF (see Table |2|, correspond to the signals TOA. In the 
present work, all selected FRT observations are binned at 460 meters per pixel (the spatial resolution of each super-pixel). 
This is done (i) to cope with the geometric deformation inaccuracies in the case of an oblique view with pointing errors, 
(ii) to minimize on poorly known topography (a CRISM pixel is smaller than MOLA resolution), (iii) to reduce local 
seasonal variations and (iv) to minimize local slopes effects. We found that binning at 460 meters is a good compromise. 



2.2 Estimation of surface bidirectional refiectance: Correction for atmospheric contribu- 
tion 

The radiative transfer in the Martian atmosphere is dominated by CO2 and H2O gases, mineral and ice aerosols which 
are an obstacle for the studies of the surface properties. Indeed the extinction of radiative fluxes, in particular the solar 
irradiance in the wavelength range of CRISM observations, is mostly due to absorption by gases and scattering by aerosols. 
Also, aerosols produce an additive signal by scattering of the solar light. As a result, a spectrum collected by the CRISM 
instrument at the Top Of-Atmosphere (TOA) is a complex signal determined by both the surface and the atmospheric 
components (gases and aerosols). The atmospheric correction chain proposed for CRISM observations is composed of: (i) 
the retrieval of Aerosol Optical Thickness (AOT), (ii) the correction for gases, and (iii) the correction for aerosols resulting 
in the estimation of the surface reflectance. The last step is carried out by the technique referred to as the MARS-ReCO 



Ceamanos et al. 20131. In the following, a short summary of the methodology is provided. 



2.2.1 Retrieval of AOT 

The AOT is defined as the aerosol optical depth along the vertical of the atmosphere layer and is related to the aerosol 
content. 

The retrieval of AOT from orbit is difficult because of the coupling of signals between the aerosols and the surface. 
Over the last decades. Emission Phase Function (EPF) sequences of Mars were obtained and were used to separate the 
atmospheric and surface contributions. Significant progresses have resulted from the pioneer works of Clancy and Lee] 



(1991 1 based on Viking orbiter InfraRed Thermal Mapper (IRTM) EPF observations and Clancy et al. (20031 from Mars 
Global Surveyor (MGS) Thermal Emission Spectrometer (TES) EPF observations ten years later. 

For our study, we decided to use Michael Wolff's AOT estimates for atmospheric correction purposes (personal com- 
munication). This parameter is available for each CRISM observation and is derived from Wolff et al. (2009)'s work. This 
method is based on the analysis of CRISM EPF sequences combined with information provided by "ground truth" results 
at both MER landing sites which allow to isolate the single scattering albedo. This method carries out a minimization of 
the mean square error between measured and predicted TOA radiance based on the previously estimated aerosol single 
scattering albedo and scattering phase function. 

Some assumptions regarding the aerosol properties (i.e., phase function, mixing ratio, column optical depth, etc) 
and the surface properties (i.e., phase function) must have been accounted for to separate the atmospheric and surface 
contributions. Concerning the surface properties, this method assumes a non-Lambertian surface to estimate the AOT 
by using a set of surface photometric parameters that appears to describe the surface phase function adequately for 
both MER landing sites ( Johnson et al. 2006a|b I . This assumption is qualitative but reasonable for several reasons 
enumerated by | Wolff et al.\ ( |2009[ ). For both MER landing sites it seems boin \Wolff et al.\ ( [2009] 's work that the AOT 
retrievals are overall consistent with optical depths returned by the Pancam instrument (available via PDS). Consequently, 
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the Wolff's AOT estimates are suitable for our study at both MER landing sites. Concerning the aerosol properties, some 
uncertainties may exist especially for the aerosol scattering phase function which is related to the aerosol particle size and 
shape. First, this method assumes a mean particle aerosol but variations are observed as a function of solar longitude, 

2006). Second the mean aerosol particle size is derived 



spatial location (latitude, longitude and altitude) ( Woljf et al. 



from CRISM observations acquired during planet-encircling dust event in 2007. During a large dust event, aerosol particle 
size are larger than those found under clear atmospheric conditions. Third, the hypothesis of the particle shape used in 
Wolff et al. (2009)'s study contain a wrong backscattering part in the scattering phase function ( Wolff et al.\ 20101). Thus 



all points previously enumerated, may bias the AOT estimation and must be taken into account during the analysis of 
the robustness of results (the surface bidirectional reflectance and the surface photometric parameters). However, the 
performance of MARS-ReCO is sensitive to the accuracy of the AOT estimate ( Ceamanos et al\ 2013| ) (see 2.2.21. 
Note that the AOT is calculated at 1 where the absorption of gases is almost null. 



2.2.2 Correction for gases and aerosols 

In order to test the performances of the methodology presented in this article and its companion paper, the present study 
is only conducted at a single wavelength. We choose to work exclusively at 750 nm where the contribution of gases is low 
and thus the retrieval of photometric properties is likely to be sufficiently accurate. Furthermore, photometric properties 
retrieved from in situ measurements taken by Pancam are available at this wavelength and our photometric properties 
can be validated. Note, however, that the presented methodology can be applied to any CRISM wavelength provided the 
contribution of gases is corrected previously. 

MARS-ReCO is devised to compensate for mineral aerosol effects considering the anisotropic scattering properties of 
the surface and the aerosols. This method is suitable for any CRISM multi-angle observation within some atmospheric and 
geometrical constraints (A0T<2, incidence angle 6'o<60°, phase angle range gmax — 5mm>40°). MARS-ReCO is based on 
a coupled surface-atmosphere radiative transfer formulation using a kernel-driven scattering model for the surface and a 



Green's function to model the diffuse response of the atmosphere (please refer to Ceamanos et al. (2013) for more detail). 
The AOT of each observation is an input of MARS-ReCO. Table [l] presents the AOT^ater and AOT mineral for each 
selected CRISM observation. We can note that AOT^ater is negligible in front of the AOT mineral and consequently, the 
photometric effects from aerosol water ice can be considered negligible in this study. 

The uncertainties pertaining to the AOT estimates (personal communication of Michael Wolff, ( Wolff et al. 2009| )) 



and to the TOA measurements by CRISM are integrated and propagated in the estimation of the surface bidirectional 
reflectance in BRF units (cf . Table [2]) . 

Besides the retrieval of the surface bidirectional reflectance, MARS-ReCO also provides an indicator of the quality 
of the estimated solution in a standard deviation sense, noted by parameter Up (in BRF units) and c omputed as — 



^ t'" (^pp)i where tr (Cpp) is the a posterior covariance matrix and Ng the available geometries 



Ceamanos et al. 



2013). In systematic test presented in the companion article, the parameter Up has proved to be highly correlated with 
the bidirectional reflectance error from MARS-ReCO and thus provides us with reliable information on the accuracy of 
the estimated surface bidirectional reflectance. 



2.3 Estimation of surface photometric properties: Bayesian inversion based on a Hapke's 
photometric model 

2.3.1 Defining regions of interest and selection 

We now discuss the improvement that results when combining different CRISM targeted observations based on their 
spatial coherence and thus a better sampling of the surface bidirectional reflectance and thus a maximization of the phase 
angle range. We defined in Subsection |2.1.3| a super-pixel as the combination of all reflectance value corresponding to a 
same location unit coming from an individual CRISM sequence or SPC cube. The combination of each super-pixel from 
each selected CRISM targeted observation is performed when their central coordinates (latitude and longitude) differ less 
than a half super-pixel size (460/2=230 meters). This is done to ensure maximum overlap. We choose same combined 
super-pixels called regions of interest (ROI) in following for the photometric study using the several criteria. First, the 
different ROIs must be located close to the MER Spirit and Opportunity rover's path, specifically to the location of 
spectrophotometry measurements by Pancam and in the same geological unit (i.e. presenting same materials). Second, 
the local topography makes the photometric study more challenging when it is poorly known because it controls to a large 
extent the incidence, emergence and azimuth local angles. Besides in the case of an oblique illumination (i.e., up to 70°), 
shadows decrease the signal/noise ratio. In this study, ROIs were therefore selected only in fiat areas. Third, ROIs are 
chosen to have the richest angular configurations and the best angular sampling in terms of phase angle range in order 
to constrain as much as possible the photometric properties. Figure |2] presents the selected ROIs for this photometric 
study at Gusev Crater and Meridiani Planum. Four ROIs are selected for Gusev Crater study while only one is chosen 
for Meridiani Planum study. The limited number of selected ROIs is explained by the fact that few ROIs in both cases 
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Figure 2: a. Image from Context camera (CTX) of the Spirit landing site at Gusev Crater with the rover path across 
the plain up to the Columbia hill in purple. The footprint of each selected CRISM observations (only the nadir image) 
is represented here (blue: FRT3192, green: FRT8CE1 and red: FRTCDA5). Full black stars point to the locations of 
the photometric measurements in the Gusev plains (1. Landing Site, 2. Bonneville Rim, 3. NW of Missoula) taken by 
instrument onboard Spirit {Johnson et al. 2006a). Full black squares represent the four ROIs that have been selected 
for our photometric study (ROI from I to IV). b. Image from CTX camera of the Opportunity landing site at Meridiani 
Planum with the rover path in purple. The footprint of each selected CRISM observations (only the nadir image) is 
represented here (blue: FRT95B8, green: FRT334D and red: FRTB6B5). Full black stars represent the locations of 
photometric measurements (1. South of Voyager, 2. Purgatory region) taken by Pancam instrument onboard Opportunity 
( Johnson et al.\ 2006b I. Full black squares represent the selected ROI for our photometric study (ROI I). 



respect the combination of criteria previously mentioned. In order to improve the number of ROIs, an improved pointing 
of each CRISM targeted observation could be envisaged. 



2.3.2 Direct surface model: Hapke's photometric model 

Models describing the photometry of discrete granular media have been proposed to express the surface bidirectional 



reflectance using semi-empirical a 
and numerical approaches (e.g.. 


analytical approaches (e.g., 


Hapke 


1981a 


Hapke' 


'1981b'; 'Hapke 


1984 


1986 1993 20021 


Cheng and Domingue\ 2000 


Doute and Schmitt^ 


1998; Mishchenko et al. 


19991. 



Hapke's modeling ( Hapke] 1993 1 is widely used in the planetary community due to the simplicity of its expression 
and due to the fact that it introduces photometric parameters claimed to have physical relevance. Previous Martian 



photometric studies have been conducted based on orbital HRSC measurements ( Jehl et al. |2008) and based on the 
analysis of in situ MER measurements ( Johnson et al. . 2006a|b I using the Hapke's 1993 version model ( Hapke 1993 1 . For 
this study, since we use the in situ investigation as "ground truth" in order to validate the MARS-ReCO approach, we 
have to use Hapke's 1993 version model. 



More recent Hapke models { Hapke\ 2002[ 2008) are available improving the original formulation {Hapke 1981a 19931 
and can be used for future photometric studies. First, the model (Hapke 2002) includes: (1) a more accurate analytic 
approximation for isotropic scatterers, (2) a better estimation of the bidirectional reflectance when the scatterers are 
anisotropic, and (3) the incorporation of coherent backscattering. Second, the model ( Hapke\ 2008) overcomes the limita- 
tions of the original model in order to predict porosity dependence of the bidirectional reflectance, in case of a particulate 
medium such as a planetary regolith. For that purpose the treatment of light propagating through the particle spacing. 
Following Johnson et al. ( 2006a|b l's works and for the sake of the comparison coherence we use the expression of Hapke 



{ 1993 1 as follows: 



{[1 + B{g)] Pig) + H{poe)HOie) - 1} 5(0o, 0, g), 



At: (/ioe + Me) 

• Geometry do, 9, and g: The incidence, emergence and phase angles, respectively. 



(1) 



• Single scattering albedo w: Factor w (0 < w < 1) depends on wavelength and represents the fraction of scattered 
light to incident light by a single particle ( Chandrasekhar 1960). 
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• Particle scattering phase function P{g)'. Function P{g) characterizes the angular distribution of an average particle. 
The empirical 2-term Henyey-Greenstein function (hereafter referred to as HG2) is used commonly for studying 
planetary surfaces (e.g., Cord et al] 2003 Hartman and Domingue 1998, Jehl et aL\ 2008 Johnson et al. 2006a|b 



Souchon et al. 2011). It reads: 



P{g) = {l-c) 

(1 + 26 cos (5) + &2)3/2 



(1-26 cos (.g) + 62) 



3/2 



(2) 



where the asymmetric parameter 6 (0 < 6 < 1) characterizes the anisotropy of the scattering lobe (from 6 = 0, which 
corresponds to the isotropic case, to 6 = 1, which corresponds to a particle which diffuses light in a single direction). 
The backscattering fraction c (0 < c < 1) characterizes the main direction of the diffusion (c < 0.5 corresponding to 
forward scattering and c > 0.5 corresponding to backward scattering). 



Multiple scattering function H{x): The exact values of the H function for isotropic scatterers was given by Chan 



drasekhar ( 


19601 


( Hapke 


19931 as 



The H93 differs by a rel ative error on H lower than 1% and a relative error on the bidirectional reflectance of a 
regolith lower than 2% [Cheng and Domingue 20001. Defining y — {1 — w^^^, 
becomes: 



the multiple scattering function 



H{x) = 1 - [1 - y] : 



1-y 
1 + y 



2 \l + y I 1 + y 



In 



1 + a; 

X 



(3) 



• Shadow Hiding Opposition Effect (SHOE) and Coherent Backscatter Opposition Effect (CBOE): B{g) is a function 
designed to model the sharp increase of brightness around the zero phase angle often observed in the case of 
particulate media, the so-called opposition effect. The B{g) function is given by Hapke\ ( "l993[ ) as follows: 



Big) = 



Bn 



l + itan(§) 



(4) 



Parameters h and Bq are respectively the angular width and the amplitude of the opposition effect. Factor h (ranging 
from to 1) is physically related to compaction and particle size distribution and Bq (ranging from to 1) is an 
empirical parameter which is related to the particle transparency [Hapke 1993). It is important to mention that the 
Ross-Thick Li-Sparse (RTLS) model for the surface reflectivity employed by MARS-ReCO ( Ceamanos et al. [2013 1 
is able to describe a backscattering lobe by means of its geometrical kernel. Nevertheless its angular width is more 
characteristic of photometric effects linked with shadows cast by macroscopic roughness than those occurring at the 



grain scale such as the SHOE and the CBOE [Lucht et al. 2000). Furthermore, the radiative transfer algorithm 



that is used to calculate the atmospheric quantities at the core of the MARS-ReCO procedure cannot propagate 
properly through the atmosphere the narrow backscattering lobes of the SHOE and CBOE. Finally, CRISM orbital 
measurements never reach the small phase angle domain (< 5°) where the previous phenomena are expressed. 
Consequently, MARS-ReCO is never in the position of retrieving accurate values for Bq and h. 

• Macroscopic roughness factor S: We note that planetary regoliths present roughness driven by grain clusters to the 
pixel scale. In the Hapke's surface model this phenomena is described by a Gaussian distribution of slopes at a 
single spatial scale under the pixel size which is not explicitly given. The mean slope angle 9 is the only required 
parameter ( Hapke\ 1993). Surface roughness involves several radiative phenomena: (i) multiple reflection of light 
between facets, (ii) shadows depending on the geometry, (iii) bias on the incidence and emergence angles, and (iv) 
increase of the multiple scattering component. In order to quantify their influence on the bidirectional reflectance, 
Hapke's model introduces a simple multiplicative factor S depending on noe and /ie whose expressions are given in 



Hapke (1993). 



2.3.3 Bayesian inversion of the surface model 

The "inverse problem" consists in estimating the model parameters that best explain the observations. Unfortunately, 
inverse problems do not have a unique solution if the direct model is non-linear, as does the Hapke's model. | Tarantola\ 
and Valette [ 1982 1 proposed to solve inverse problems in a general non-linear case based on the concept of the state of 



information which is characterized by a Probability Density Function (PDF). The PDF is defined over both the parameter 
space and the observed space. The formalism of a PDF is used to define the initial state of information (i.e., a priori 
knowledge on the parameters, the uncertainties on the observation and on the model). To infer the solution, the Bayes' 
theorem is applied. Key points concerning the bayesian inversion concept and framework assumptions are presented in 
the following: 
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• Data, model parameters and theoretical relationship: The direct model consists of computing the simulated data d, 
from model parameters m: 



d = F{m) 



(5) 



• Prior information on the model: The prior information on model parameters in the parameter space (M) 
is independent with the data and corresponds to the state of null information. For the Hapke model parameters 
uj, b, c, 0, Bq and h, we consider a uniform PDF, different from zero on an interval that insures their physical 
relevance (between and 1 for lj, b, c, and h and between 0° and 90° for 0). Outside the intervals, the PDF 
is null, avoiding unphysical solution to appear. As discussed in Subsection |2.3.2| in real planetary situations, the 
SHOE and the CBOE phenomena are only expressed for phase angles 5 < 5° out of the range encompassed by 
typical CRISM observations {g > 30°). Thus a priori, neglecting both phenomena should not influence the retrieval 
of the other parameters . However following Souchon et al. (20111, the model can still be profitably inverted on 
the data by keeping parameters Bq and h since they will compensate for discrepancies between the model and the 
measurements in some situations. Consequently, in the present work, we first chose to invert the parameters Bq and 
h to systematically control if they are constrained or not by the CRISM data set. However in a second phase we 
also tested the inversion by setting Bq and h to zero and no change was observed on the determination of the other 
parameters (w, 6, b, c). 

• Prior information on the data: The prior information on data pdid) in the observation space (D) is assumed to 
be a Gaussian PDF according to the MARS-ReCO formalism and retrieval strategy. Note that the error tr on a 
CRISM measurement at one geometry is assumed to be independent on the state of the surface and on the other 
geometries (i.e., the PDF has a diagonal covariance matrix with elements (T^,...,cr^ , where Ng is the number 
of available geometries, up to 11). At 750 nm, the signal to noise ratio was estimated before launch to be equal 
to 450 (Murchie et al. 2007[ ) but due to additional artifacts such as spikes, and calibration issues {Seelos et al. 
2011 1 we evaluated the uncertainty a of the reflectance measurement (p) at each geometry j to be o f the order 

is the CRISM dataset at the j*'* angular configuration ( | Ceamanos 



of cr, = (1/50) X pj where j=l. 



.Ng and 



et al. |2013[ ). Moreover, MARS-ReCO takes into account the uncertainty of the AOT input. Figure [s] presents a 



typical TOA photometric curve collected by the CRISM instrument (green plus) and the bidirectional refiectance 
curve produced by the MARS-ReCO algorithm (red crosses). For each geometry the bidirectional reflectance value 
is accompanied by its la uncertainty. Those means and root mean square errors are used to build pdid) that serves 
as an input PDF of the Bayesian inversion. 

• Posterior probability density function and resolution of inverse problems: Inversion problems correspond to the 
particular case where information from the data space (D) is translated into the model space (M). The posterior 
PDF in the model space (jM{fn) as defined by Tarantola and Valette (19821 reads 



<JMim) = kpM(m) L(m), 
where fc is a constant and L{m) is the likelihood function 

PD{d)9id\m) 



L{m) = / dd 
Jd 



paid) 



(6) 



(7) 



where 0{d \ m) is the theoretical relationship of the PDF for d given m, and poid) is null information PDF for the 
data d. The solution of the general inverse problem is given by the PDF aM (m) from which any information on the 
model parameter can be obtained such as mean values, uncertainty bars, etc . Please refer to \Tarantola and Valette\ 



( 1982 ) for more detail. 



• Sampling of solutions to inverse problems: In our case, the relationship between model parameters and observed 
data through Hapke's modeling is non-linear. While it is not possible to describe the posterior PDF analytically, 
it can be sampled using a Monte Carlo approach using a Markov chain (Mosegaard and Tarantola 19951. After 



a sufficient number of steps, the state of the chain corresponds to the desired distribution. According our tests, 
the best trade-off between computation time and accuracy is a burn-in phase (phase in which the Markov chain 
approaches a stationary state after a certain number of runs) of 500 steps. The next 500 iterations are used to 
estimate the posterior PDF allowing the determination of the mean and standard deviation of each parameter. Note 
that a posteriori PDF of a retrieved parameter is not necessarily a gaussian distribution but can be a square function 
or multi-modal as seen in Figures |4] and |6] To describe the results of each parameter, we choose to compute the 
mean of the posterior PDF. To describe the uncertainties, we choose to compute the standard deviation. We warn 
the reader because sometimes these estimators can be inappropriate to describe the PDF. In the following graphs, 
2(7 error bars are plotted to describe more accurately the PDFs. 
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Figure 3: Reflectance values (in BRF units) corresponding to the photometric curve of the ROI I from FRT3192 observation 
composed of ten angular conflgurations at the TOA (green plus), and at the surface after the correction for atmospheric 
contribution carried out by MARS-ReCO (red crosses). In both cases the la uncertainty is expressed as a function of 
angular configuration (the eleven geometric views of one CRISM EPF scans from 1st to 5th corresponds to the inbound 
direction, 6th is the central scan and from 7th to 11th are related to the outbound direction). The latter data are used 
in the Bayesian inversion. The blue solid and dashed lines represent the Hapke's best fit and their la uncertainties, 
respectively. 



• Root Mean Square residual RMS: In order to estimate the difference between the fit and the observed bidirectional 
refiectance, the Root Mean Square residual (noted RMS) is given for each Bayesian inversion as follows: 



RMS 



X^ATg {Pobs Pmod) 



where N is the available geometric configurations, Pobs the CRISM bidirectional refiectance corrected for atmosphere 
and Pmod the modeled bidirectional refiectance taken as the mean of the 500 iterations used to estimate the posterior 
PDF. 

Non-uniformity criterion k: Photometric parameters m are constrained if their marginal posterior PDF differs from 
the prior state of information (i.e., a null information taken as a uniform distribution, in our case). In order to 
distinguish if a given parameter has a solution we perform a statistical test leading to a non-uniformity criterion k 
(see Appendix A). For k > 0.5, the marginal posterior PDF is considered to be non-uniform and thus we consider 
that the mean and standard deviation of the PDF satisfactorily describe the solution(s). 



3 Analysis of retrieved photometric parameters 

As mentioned in Subsection |2.1.2| we selected three CRISM targeted observations from Gusev Crater and from Meridiani 
Planum (cf. Table [ij. All observations were individually corrected for aerosols using MARS-ReCO, knowing the respective 
AOT values (cf. Table [ij. Based on the corrected bidirectional refiectance, we use the proposed methodology described in 
Subsection |2 . 3 . 3 1 to estimate the photometric parameters of the materials encompassed by four ROIs in the case of Gusev 
Crater (ROI I to ROI IV) and by one ROI in the case of Meridiani Planum (ROI I). 

For each ROI, we determine the surface photometric parameters at 750 nm from (i) a single CRISM FRT observation 
or (ii) a combination of all selected FRT observations (i.e., FRT3192, FRT8CE1 and FRTCDA5 in the case of Gusev 
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Table 3: Retrieved Hapke's parameters {lo, b, c, 6 in degree, Bq, h) and their standard deviation in parentheses (the 
underconstrained parameters are indicated by "(-)") corresponding to the ROI I and II from CRISM measurement at 750 
nm using our Bayesian inversion as well as the non-uniform criterion k (second line of each entry) after MARS-ReCO. This 
analysis is done for each CRISM FRT observation (i.e., FRT3192, FRT8CE1 and FRTCDA5) or for the combination of the 
three observations after correction for the atmospheric effects by MARS-ReCO. The standard deviation a-p (in BRF units) 
is a parameter provided by MARS-ReCO estimating the quality of the retrieved surface bidirectional reflectance. The 



absolute quadratic residual of the flt RM S (in BRF units) evaluates the deviation of the flt from observed bidirectional 
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0.33±0.Q4 


0.01 
0.01 


0.68 (0.06) 
1.00 


0.17 (0.20) 
1.92 


0.62 (0.20) 
0.96 


11.62 (3.98) 
1.00 


0.52 (-) 
0.24 


0.52 (-) 
0.12 
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8CE1 


-14.606 
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0.98±0.15 


0.03 
0.01 


0.77 (0.07) 
1.00 


0.30 (0.28) 
1.47 


0.42 (0.27) 
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17.23 (9.64) 
0.98 


0.47 (-) 
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0.50 (-) 
0.12 
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0.01 
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0.51 (-) 
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0.01 
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1.00 
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0.23 


0.52 (-) 
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0.01 


0.72 (0.05) 
1.00 


0.27 (0.15) 
1.18 


0.56 (0.16) 
0.99 


15.62 (2.43) 
1.00 


0.49 (-) 
0.05 


0.52 (-) 
0.23 
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Crater, and FRT95B8, FRT334D and FRTB6B5 in the case of Meridiani Planum). Remember that the chosen targeted 
observations are complementary in terms of phase angle range (cf. Table [l}. 

For each parameter of the photometric Hapke's model (i.e., uj, b, c, 9, Bq and h) we determine its mean value and 
standard deviation by running the proposed Bayesian inversion procedure. The non-uniform criterion noted k is conjointly 
computed and detailed in Tables |3] and |4] for Gusev Crater and Table [5] for Meridiani Planum. In the following, we have 
two criteria that help us estimate the existence and the quality of a solution. On the one hand, the non-uniform criteria 
k allows us to reject posterior PDF's that do not carry a solution. On the other hand, the standard deviation allows us 
to quantify to which extent a given parameter is constrained by the current solution. 

3.1 Results on Gusev Crater 

The quality of the surface bidirectional reflectance estimated by MARS-ReCO is given by the standard deviation dp. This 
quality parameter is available for each ROI of each CRISM targeted observation of the present study (cf . Tables |3] and |4]) . 
The highest a-p value is observed for FRT8CE1 which can be explained by a wrong AOT estimation (highest AOT value) 
in this case (i.e., FRT8CE1: ap = 0.03-0.04 for AOT = 0.98). The positive correlation of this uncertainty and the AOT 



is also observed in the sensitivity study led by Ceamanos et at. (20131. Indeed the error computed for synthetic reference 
data mimicking the photometric properties of the planet Mars increases with AOT (i.e., for AOT = 1, dp ~ 0.05 and error 
is ^10% while for AOT = 1.5, dp ^ 0.10 and error is ~20%, see Ceamanos et al. (20131 for more detail). 



Tables [3] and [4] present results regarding the Hapke's model parameters for the different ROIs of Gusev Crater. 
The goodness of flt is estimated through the absolute quadratic residual RMS value (cf . Tables |3] and |4]) . For all 
Bayesian inversions, the estimates are less than 0.02 which mean that the inversions are deemed satisfactory. 

Figure |4] represents the PDF of each parameter considering the inversion of a single FRT observation as well as the 
inversion of the three combined observations. According to results, a solution exists for parameters uj, b, c and 9 (because 
the PDF is non-uniform) whereas no solution exists for parameters Bq and h (uniform PDF). Several conclusions can be 
drawn when using a single targeted observation: 

• Solutions exist for parameter lu in all cases (fc ~ 1). The standard deviation shows that the single scattering albedo 
is the most constrained parameter (i.e., 0.06 < a < 0.09). Examples of posteriori PDFs are plotted in Figures |4]and 

m 
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Table 4: Same as Table [3l but for the ROI III and ROI IV 
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• Although solutions exist for parameter h in all cases (fc » 0.5), the standard deviation shows that it is poorly 
constrained (i.e., 0.25 < a < 0.32) except for the ROI I and ROI IV using the CRISM observation FRT3192. This 
discrepancy can be explained by the higher number of available geometric configurations, respectively 10 and 9. 
Examples of posteriori PDFs are plotted in Figures |4] and [TO] 

• We find meaningful values for parameter c only when we use FRT3192 or FRTCDA5 (fc > 0.5). The standard 
deviation is relatively low (i.e., 0.20 < a < 0.23) in this case. In the case of FRT8CE1, no solution is found for 
parameter c {k < 0.5) for ROI II, III and IV. For the ROI I, however, a solution exists but it is poorly constrained 
(i.e., a = 0.27). Examples of posteriori PDFs are plotted in Figures |4] and 10 



• While parameter 9 has a solution in all cases {k ^ 1), we distinguish two types of results: (i) using FRT3192 or 
FRTCDA5, we note that the standard deviation is relatively low (i.e., 3.98 < a < 5.81), and (ii) using FRT8CE1, 
the standard deviation is relatively high (i.e., 9.64 < a < 20.40) .Examples of posteriori PDFs are plotted in Figures 
Handllll 



• No solutions are found for parameters Bq and h in any case {k « 0.5). As mentioned in Subsection 2.3.2 



nor 



the phase domain covered by the CRISM observations nor the capabilities of MARS-ReCO allow to constrain the 



opposition effect {Ceamanos et at. 20131. Consequently, accurate values for Bq and h cannot be retrieved and no 



physical interpretation can be done in this case. Examples of posteriori PDFs are plotted in Figure |4] 

Regarding the processing of a single CRISM observation, we note that parameters c and 9 are, respectively, non- 
constrained and poorly constrained only when treating data from FRT8CE1. The reason to explain such a difference may 
be that the available maximum phase angle range is less than 90° for FRT8CE1. Indeed, \Helfenstein\ (|1988|) underhned 



the necessity to have observations which extend from small phase angles out to phase angles above 90° for an accurate 
determination of the photometric roughness. By contrast, the phase angle range expands by more than 100° for the 
other observations (cf . Table [T]) . We note that parameter b is poorly constrained when treating data from all available 
observations. This result may be explained again by the available phase angle range (cf. Table [T]). In conclusion, the 
presented study clearly demonstrates that the existence and quality of a solution for parameters b, c and 9 is dependent 
on the available phase angle range. The bidirectional reflectance curve from a single CRISM observation does not contain 
enough phase angle information. 

By contrast, the combination of three targeted observations provides improved constrains on all Hapke's parameters 
except for Bq and h. Indeed, the standard deviation of each estimated parameter is lower than those obtained when using 
only a single targeted observation. We note that for ROI IV, the Hapke's parameters are less constrained than for the 
other ROIs, especially for parameter b for which no solutions were found. This result can be explained by the limitation in 
the bidirectional sampling. In fact, we can note in Figure [5] which represents the north projection of geometric conditions 
for the four selected ROIs that the ROI IV misses a near nadir geometry from the observation FRT8CE1. The latter 
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Figure 4: The PDF of each of the six physical parameters retrieved for ROI I. Two cases are considered: the inversion of a 
single CRISM FRT observation (i.e., FRT3192 in blue line) and the inversion of three combined CRISM FRT observations 
(in red line). 



FRT 


Lat. 


Lon. 


AOT 






b 


c 


e 


Bo 


h 


nh 


9 




(°N) 


(°E) 




RMS 


k 


k 


k 


k 


k 


k 






95B8 


-2.000 


-5.484 


0.56±0.09 


0.02 
0.01 


0.72 (0.09) 
0.99 


0.34 (0.30) 
1.40 


0.42 (-) 
0.44 


17.32( 10.33) 
0.98 


0.48 (-) 
0.12 


0.50 (-) 
0.16 


5 


~41-86 


334D 


-2.001 


-5.482 


0.35±0.04 


0.01 
0.01 


0.67 (0.09) 
0.99 


0.38 (0.32) 
1.16 


0.49 (0.25) 
0.61 


16.65 (7.82) 
1.00 


0.49 (-) 
0.18 


0.48 (-) 
0.16 


5 


~50-98 


B6B5 


-2.004 


-5.482 


0.35±0.04 


0.01 
0.01 


0.65 (0.07) 
1.00 


0.19 (0.20) 
2.06 


0.52 (-) 
0.41 


18.09 (5.82) 
1.00 


0.50 (-) 
0.12 


0.50 (-) 
0.14 


6 


~41-106 


ROI I - 3 FRTs 








0.02 


0.68 (0.08) 
0.99 


0.36 (0.26) 
0.82 


0.41 (0.21) 
0.95 


18.37 (5.20) 
1.00 


0.51 (-) 
0.15 


0.52 (-) 
0.24 


16 


~41-106 



shows a more different incidence angles (nearly 40°) than the two other CRISM FRT observations (nearly 60°) (see also 
Table [T]) thus explaining that the parameters are less constrained for ROI IV than the other ROI. 

To conclude, the quality of these results show the benefits of combining several targeted observations in constraining 
photometric parameters. 

3.2 Results on Meridiani Planum 

Similar to the photometric study on the Gusev Crater, the standard deviation o-p determined by MARS-ReCO is given 
for each CRISM observation and the ROI I used for the present study (cf. Table [S]). Note that dp values are acceptable 
(0.01 < CTp < 0.02) meaning that estimated surface bidirectional reflectance are accurate. 

Table |5] presents the results obtained on Meridiani Planum. Similar to Gusev Crater study. Figure |6] shows that a 
solution exists for the parameters u, b, c and 6 (i.e., non-uniform PDF) when a single CRISM FRT and in the case of 
three combined observations, except for parameters Bq and h. 

The goodness of fit is estimated through the absolute quadratic residual RMSabs value (cf. Table |5|. For all bayesian 
inversions, the estimates are less than 0.02 which mean that the inversions are deemed satisfactory. 

We note that for the parameter lo, two maximum are visible at nearly 0.6 and 0.8. In order to understand the origin of 
the bimodal distribution for the parameter uj in this case, the refiectance of typical Martian material was generated by using 
realistic photometric properties determined by the Pancam instrument aboard the MER Opportunity site in Meridiani 
Planum (i) at the same geometric configurations as the FRT95B8 observations, and (ii) at the same combined geometric 
configurations when merging the three selected observations. In both cases, the posterior PDF for the parameter lu shows 
a bimodal distribution. If the geometric sampling is broader (i.e., with varied incidence, emergence and phase angles), the 



12 




Figure 5: North projection of geometric conditions (stars: incidence rays and crosses: emergence directions) of the four 
selected ROIs for the Gusev Crater study (blue: FRT3192, green: FRT8CE1 and red: FRTCDA5). 



posterior PDF of the parameter u becomes a gaussian with a single peak. The presence of two possible solutions is the 
consequence of the limitation of a sufficient geometric diversity in our selection of CRISM observations for the Meridiani 
Planum study to constrain the parameter cj, which is otherwise the best-constrained parameter in photometric modeling. 
Several conclusions can be drawn when using a single CRISM targeted observation: 

• Solutions exist for parameter uj in all cases (/c ~ 1). The standard deviation shows that the single scattering albedo 
is the most constrained parameter (i.e., 0.07 < a < 0.15). Examples of posteriori PDFs are plotted in Figures [6] and 

[HI 

• We find meaningful values for parameter h only when for FRT95B8 or FRT334D or FRTB6B5 (k > 0.5). Albeit 
the standard deviation is relatively high for FRT95B8 or FRT334D (i.e., 0.30 < a < 0.32), the standard deviation 
becomes relatively low for FRTB6B5 (i.e., a = 0.20). Examples of posteriori PDFs are plotted in Figures |6]and[l^ 



Solutions only exist for parameter c when using FRT334D (fc > 0.5). However, the standard deviation is relatively 



14 



high (cr = 0.25). Examples of posteriori PDFs are plotted in Figures [6] and 

• While parameter 6 has a solution in all cases {k ^ 1), we distinguish two types of results: (i) using FRT334D or 
FRTB6B5, we note that the standard deviation is relatively low (i.e., 5.82 < a <7.82) and (ii) the standard deviation 
is relatively high (i.e., a = 10.33) when using FRT95B8. Examples of posteriori PDFs are plotted in Figures [6] and 

m 

• Similar to the Gusev study, no solution is found for parameters BO and h in any case {k » 0.5). Examples of 
posteriori PDFs are plotted in Figure [6j 

Dealing with single observations, we note that parameter 6 is poorly constrained when treating data from FRT95B8 
and FRT334D whereas it is highly constrained when using FRTB6B5 which can be explained by a maximum phase 
angle below 90° in case of FRT95B8 and FRT334D (cf. Table [l]). \Helfenstein\ ^988) underHned the necessity to have 



observations which extend from small phase angles out to phase angles above 90° for an accurate determination of the 
photometric roughness. Note that parameter b and c are non-constrained or poorly constrained even through solutions 
exist in all cases. This outcome can be explained by a worse quality of the Meridiani Planum bidirectional reflectance 
sampling. Indeed, we note in Figure [7] which represents the north projection of geometric conditions of each selected 
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Figure 6: The PDF of each of the six physical parameters retrieved for ROI I of Meridiani Planum. Two cases are 
considered: the inversion of a single CRISM FRT observation (i.e., FRT95B8 in blue line) and the inversion of three 
combined CRISM FRT observations (in red line). 

CRISM FRT observations that the three selected CRISM FRT (i) miss near nadir geometries (close to emergence equal 
0) and (ii) present lower number of available angular configurations compared to Gusev Crater work (lower than 6). 

We then improve the bidirectional reflectance sampling by combining the three selected targeted observations. We 
observe that parameters uj, b, 9 become significantly more constrained. Indeed, the standard deviation of each estimated 
parameter is lower than those obtained when using only a single observation. However, the parameters are less constrained 
than those obtained for the Gusev Crater study, which can be explained by the lack of near nadir geometry for this case 
(cf. Figure [7]). 

Again, results underline the benefits of combining several observations. 



4 Validation 

This section focuses on the validation of the estimated photometric parameters by comparing them to previous photometric 
studies based on experimental, in situ and orbital photometric studies. As the PDF of each parameter is not necessary 
a gaussian, the estimated mean and the standard deviation are not always representative of the entire distribution. 




Figure 7: North projection of geometric conditions (stars: incidence directions and crosses: emergence directions) of the 
selected ROI for the Meridiani Planum study (blue: FRT95B8, green: FRT334D and red: FRTB6B5). 
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Consequently for each parameter, we use the PDF built from the last 500 iterations of the Markov chain process (see 
2.3.1 1 instead of the previous statistical estimators for the comparison to previous studies. 



4.1 Validation of results from Gusev crater 

4.1.1 Comparison to experimental measurements on artificial and natural samples 



McGuire and Hapke { 1995 1 studied the scattering properties of different isolated artificial particles which have different 
structure types (sphere/rough particles, clear/irregular shape particle, with/without internal scatterers, ...). Their study 
showed that the Henyey-Greenstein function with two parameters, HG2 (backscattering fraction c and asymmetric pa- 
rameter b) provided the best description of their laboratory bidirectional reflectance measurements. In short, their study 
shows that smooth clear spheres exhibit greater forward scattering (low values of c) and narrower scattering lobes (high 
values of b) whereas particles characterized by their roughness or internal scatterers exhibit greater backward scattering 
(high values of c) and broader scattering lobe (low values of b). In a graph mapping the b and c parameter space, the results 
exhibit a "L-shape" from particles with high density of internal scatterers to smooth, clear, spherical particles. For their 



study, McGuire and Hapke { 1995 ) used centimeter sized artificial particles which are larger than the light wavelength and 



far than typical constituents of the planetary regoliths. In order to examine the impact of particle size on the scattering 



phase function, Hartman and Domingue (1998) observed that there are no significant variations of the latter in McGuire 
and Hapke ( 1995 ) measurements when the particle size is similar to typical planetary regoliths particles. Hartman and 



Domingue (1998) concluded that McGuire and Hapke (19951 results could be considered to be representative of their 



respective particle structure types independent of particle size. However, recent experimental works have questioned the 
initial interpretation of the Hapke's parameters. Indeed, the Hapke's parameters seem to be more sensitive to the orga- 
nization of the surface material (effects of close packing, micro-roughness) than to the optical properties depending on 
individual particles ^Cord et al. [2003] Shepard and Helfenstein\ 2007). Similarly, Souchon et al. (2011) measured for a 
comprehensive set of geometries the refiectance factor of natural granular surfaces composed of volcanic materials differing 
by their grain size (from the micron-scale to the millimeter-scale), shapes, surface aspect, and mineralogy (including glass 
and minerals). Thus the main novelty of Souchon et al. ( 2011[ )'s experimental study compared to McGuire and Hapke 
(1995 1 's works is the determination of the parameters b and c for planetary analogs of basaltic granular surfaces and not 



isolated artificial particles. Souchon et al.\ (2011) compared the scattering parameters retrieved by inversion of Hapke's 



model with results on artificial materials \McGuire and Hapke 1995 1 and a similar trend was found, though with some 
variations and new insights. Granular surfaces even with a moderate proportion of isolated translucent monocrystals 
and/or fresh glass exhibit strongly forward scattering properties and a new part of the L-shape domain in the b and c 
parameter space is explored. 

In Figure [s] that is given as an example of similar figures, the scattering parameters (i.e., backscattering fraction c and 
asymmetric parameter b) of the ROI I from CRISM data (cf. Table [3| are plotted along with the scattering parameters 



obtained from laboratory measurements of artificial (McGuire and Hapke 1995) and natural samples (Souchon et al 



2011). Results shows that parameters b and c retrieved from the inversion of combined CRISM FRT observations are 
consistent with the laboratory studies. The combination of three FRTs is necessary to constrain satisfactorily the 6 and 
c values with acceptable error bars. As it can be seen, surface materials of ROI I have high values of c and low values 
of 6, which indicate broad backscattering properties related to artificial materials composed of spheres with a moderate 
density of internal scatterers and close to irregular or round rough, opaque and solid natural particles (cf. Figure Isl). 



4.1.2 Comparison to in situ measurements taken by Pancam/MER-Spirit 

The Pancam instrument on-board MER-Spirit acquired several spectrophotometric observations along the rover's traverse 
paths to determine the surface physical and chemical properties of rocks and soils encountered at the Gusev Crater. 
Johnson et al. (2006a I evaluated the parameters of Hapke's photometric model (Hapke 1993) using Equations [l] [2] and 



3| from measured radiance for several identified units (i.e., "Gray" rocks, which are free of airfall-deposit dust or other 
coatings, "Red" rocks, which have coating, and "Soil" unit, corresponding to unconsolidated materials). The measured 
radiance at the ground was first corrected for diffuse sky illumination (Johnson et al. 2006a) and local surface facet 



orientations (Soderblom et al. 2004). 



The difference of spatial resolution between CRISM and Pancam instruments must be considered prior to comparison. 
Indeed, a direct comparison between both sets of photometric parameters must be handled with care as Pancam provides 
local measurements (at a centimetric scale such as soils and rocks) , whereas orbital instruments such as CRISM measure 
extended areas integrating different geologic units (at a pluri-decametric scale such as a landscape). In the latter case, 
measurements may be dominated by unconsolidated materials such as soils. As mentioned in Subsection |2.3.1| the four 
selected ROIs are chosen close to the MER-Spirit's path. The retrieved photometric parameters for each ROI from the 
combination of three CRISM observations were compared to those extracted from three sequences taken by Pancam in 
the Gusev plains. These measurements are located near to the four ROIs in Landing site (Sol 013), Bonneville Rim (Sol 
087-088) and Missoula (Sol 102-103) to the northwest side of Columbia Hill (cf. Figure [2| ^Johnson e^"al| |2006a[ ) . In this 
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Figure 8: Probability density map of the asymmetric parameter b (horizontal axis) versus the backscattering fraction c 
(vertical axis) solutions retrieved at 750 nm estimated from the last 500 iterations of the Bayesian inversion for FRT3192, 
FRT8CE1 and FRTCDA5 observations and their combination (ROI I). The grid is divided into 24 (vertical) x 20 (hor- 
izontal) square bins. The coloring gives the probability corresponding to each bin. Means and 2a uncertainties (red 
rhombus) are plotted too. The inversion solutions are plotted against the experimental b and c values pertaining to 
artificial particules measured by McGuire and Hapke (19951 and to natural particles measured by Souchon et al. (2011). 
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Figure 9: Mean and uncertainties of the single scattering albedo uj estimated from Pancam measurements at 753 nm 
for different geo l ogical units at Landing Site (Sol 013), Bonneville Rim (Sol 087-088) and NW of Missoula (Sol 102-103) 
{ Johnson et al.\ 2006a I compared to those estimated from CRISM measurements at 750 nm (note that the error bar 
represents 2a uncertainties) derived from 2-term HG models. The PDFs of the parameter lu estimated from the last 500 
iterations of the bayesian inversion are also represented on the left side. This is helpful when mean and uncertainties 
are not entirely representative of the PDF (such as the PDF of the parameter ui for the ROI IV which shows a bimodal 
distribution) . 



case the combination of three CRISM observations were treated using the H93 version of the Hapke's model in order to 
match the same model used in Johnson et al. (2006a I. Pancam results are summarized in Table [6] 



Table 6: Retrieved Hapke's parameters and their standard deviation (underconstrained parameters are indicated by "(+,- 
)") from the 2-term model for Sol 13, Sol 87-88 and Sol 102-103 and each of the three considered units (Gray Rock, Red 
Rock and Soil) at 753 nm (except for the parameter c of Red Rock unit for Landing Site which corresponds to the 754 



nm model results) from Pancam onboard Spirit ( Johnson et al. 2006a ) 



Site 


Unit 








e (deg 


) 




b 




c 




Number 


9 (deg.) 


Landing Site 


Bright Soil 


0.75 ( 


+0.01, 


-0.00) 


2 ( 




) 


0.243 ( 


+0.020, 


-0.017) 


0.625 (+0.012, 


-0.013) 


21 


~30-120 


(Sol 013) 


Gray Rock 


0.83 ( 


+0.01, 


-0.01) 


7 ( + 


3, 


4) 


0.931 ( 


+0.045, 


-0.044) 


0.065 (+0.058, 


-0.058) 


66 


~30-120 




Red Rock 


0.79 ( 


+0.03, 


-0.02) 


20 (- 


h3. 


-3) 


0.187 ( 


+0.026, 


-0.031) 


0.720 (^0.084, 


-0.087) 


68 


~30-120 




Soil 


0.76 ( 


+0.01, 


-0.01) 


15 (- 


h2. 


-1) 


0.262 ( 


+0.010, 


-0.010) 


0.715 (+0.029, 


- 0.032) 


51 


~30-120 


Bonneville Rim 


Gray rock 


0.72 ( 


+0.04, 


-0.04) 


23 (- 


h3. 


-4) 


0.434 ( 


+0.035, 


-0.037) 


0.359 (-0.048, 


-0.050) 


63 


~25-120 


(Sol 087-088) 


Red rock 


0.70 ( 


+0.01, 


-0.01) 


15 (- 


h3. 


-3) 


0.219 ( 


+0.017, 


-0.020) 


1.000 (+, 


-) 


72 


~25-120 




Soil 


0.66 ( 


+0.00, 


-0.00) 


7 ( + 


1, 


1) 


0.170 ( 


+0.008, 


-0.008) 


0.823 (-0.025, 


-0.024) 


54 


~25-120 


NW of Missoula 


Gray rock 


0.70 ( 


+0.02, 


-0.02) 


13 (- 


h2. 


-3) 


0.406 ( 


+0.018, 


-0.032) 


0.206 (-0.019, 


-0.025) 


225 


~0-125 


(Sol 102-103) 


Red rock 


0.83 ( 


+0.02, 


-0.02) 


19 (- 


hi. 


-2) 


0.450 ( 


+0.023, 


-0.064) 


0.255 (+0.032, 


-0.072) 


141 


~0-125 




Soil 


0.69 ( 


+0.01, 


-0.00) 


11 (- 


hi. 


-1) 


0.241 ( 


+0.011, 


-0.009) 


0.478 (^0.036, 


-0.022) 


132 


~0-125 



Single scattering albedo: Figure |9] plots means and uncertainties as well as the PDF of single scattering albedo 
ui estimated for the four selected ROIs against values extracted for each unit defined by Johnson et al. (2006a I 
(i.e., Gray rocks, Red rocks and Soil) at the Landing site, Bonneville rim and NW of Missoula. The bimodality is 



identical to Meridiani ROI I (see subsection 3.2) and is the consequence of the lack of geometric diversity in the 



CRISM dataset. The estimated values of uj from CRISM are consistent with the Pancam outputs for the Soil unit. 

• Phase function: Figure [TO] plots the asymmetric parameter b versus the backscattering fraction c solutions provided 
by the Bayesian inversion for each ROI in comparison to the photometric parameters obtained by \Johnson et al.\ 
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Figure 10: Probability density map of the asymmetric parameter b (horizontal axis) and the backscattering fraction c 
(vertical axis) solutions from the last 500 iterations of the Bayesian inversion estimated at 750 nm obtained for the four 
selected ROIs. The grid is divided into 24 (vertical) x 20 (horizontal) square bins. The coloring gives the probability 
corresponding to each bin. Means and 2cr uncertainties (red rhombus) are plotted too. The inversion solutions are plotted 
against the experimental b and c values pertaining to artificial particles measured hy \McGuire and -ffapfce| (|1995^ and to 
natural particles measured hy ^Souchon et aZ.| (2011l. The Gray Rock (circle), Red Rock (square) and Soil (triangle) units 
from Pancam sequences at Landing Site (green), Bonneville Rim (blue) and NW of Missoula (orange) estimated at 753 
nm (except for the Red Rock unit at the Landing Site which is constrained at 754 nm) ( Johnson et al.\ 2006a) are plotted 
here. 



(2006a I. In the case of Red Rock unit for Landing Site, the parameter c is determined at 754 nm instead of usual 
753 nm. Moreover, the parameter c is not constrained at both wavelengths in the case of Red Rock at Bonneville 
rim. Consequently the parameters b and c are not plotted here. The b and c pair from CRISM observations are 
most consistent with the Soil unit for all the studied Spirit sites and the Red Rock unit from the Landing Site 
which exhibit broad backscattering properties. Again, we note that for ROI IV, the Hapke's parameters are less 
constrained than for the other ROIs, especially for parameter b for which no solutions were found. 

• Macroscopic roughness: Means and uncertainties as well as the PDF of the macroscopic roughness 9 obtained for each 
selected ROI are presented in Figure 11 The parameter 6 of all selected ROIs, estimated from CRISM measurements 
is consistent with the Soil unit found at the Landing Site and at NW of Missoula (cf. Figure 11). By contrast, this 



parameter indicates a rougher surface than the soil unit found at Bonneville Rim which is consistent with \Johnson\ 



et al. (2006a I. At intercrater plains where the Landing site is situated, the greater proportion of small clasts in the 



soil (compared to the soils near Bonneville rim and Missoula area) may explain the rougher surface texture ( Ward 

" TZ|po55| . 



Two principal points from the previous comparison can be outlined. First, at the effective spatial scale achieved from 



space by CRISM for the photometric characterization (460 m), the surface behaves like the Soil units defined by Johnson 



et al. (2006a) and particularly like the Soil unit observed at the Landing Site. Second, the values of lu and (6, c) resulting 
from our analysis are approximately inside the range of variation seen in Figures |9] and [TO] for the Soil units at the three 
locations. This is not the case for the macroscopic roughness which is close to (within the error bars) the Soil unit values 
only for the Landing Site and Missoula area (cf. Figure 111. Consequently all photometric parameters retrieved following 



the methodology proposed in this article are consistent with the results arising from the analysis of Pancam measurements 
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Figure 11: Mean and uncertainties of the macroscopic roughness estimated from Pancam measurements at 753 nm 
for different geo logical units at Landing Site (Sol 013), Bonneville Rim (Sol 087-088) and NW of Missoula (Sol 102-103) 
( Johnson et al.\ |2006aJ compared to those estimated from CRISM measurements at 750 nm (note that the error bar 
represents 2a uncertainties) derived from 2-term HG models. The PDFs of the parameter 9 estimated from the last 500 
iterations of the Bayesian inversion are also represented on the left side. This is helpful when mean and uncertainties are 
not representative of the PDF. 



at Gusev provided that ROIs are associated to the proper soil unit. That means that the two independent investigations 
(the present work and the Johnson et al. ( 2006a[ )'s study) are cross- validated in the Gusev case. 



4.1.3 Comparison with photometric estimates derived from orbital measurements by HRSC/MEx 



Using HRSC imagery, which provides multi-angle data sets (up to five angular configurations by orbit) (McCord et al 



2007 ) , Jehl et al. { 2008 ) determined the regional variations of the photometric properties at the kilometer spatial scale 
across Gusev Crater and the south flank of ApoUinaris Patera using several observations in order to cover a phase angle 
range from 5° to 95°. The photometric study was carried out without any atmospheric correction but ensuring that the 



atmospheric contribution was limited by selecting HRSC observations with AOT lower than 0.9. Jehl et al. (2008) applied 
an inversion procedure developed by Cord et al. (20031 based on the H93 version of Hapke's model (equation [1| . In Jehl 
et al. (2008 1 's study, even photometric units were determined using a principal component analysis at 675 nm in which 
one of them corresponds to the Spirit landing site area. There is a robust overall first-order consistency between these 
photometric orbital estimates independently retrieved from HRSC and CRISM observations, particularly for ROIs I, II 
and III (see Table [7|. 

The principal point from the comparison between both orbital photometric results is that there is a robust overall first- 
order consistency between these photometric orbital estimates reached independently from HRSC and CRISM observations, 
particularly for ROIs I, II and III (see Table [T]). However, in detail, two differences can be noted for the parameters oj 
and c (see Table [t]). In fact, the parameters ui and c from CRISM measurements are respectively lower and higher than 
those determined from HRSC measurements. Moreover, the parameters uj and c estimated from CRISM data set are more 
consistent with in situ measurements. These differences are explained by the use of an aerosol correction in the CRISM 
data whereas no correction was made in HRSC data. As a consequence, the contribution of bright and forward scattering 
aerosols in the HRSC measurements appeared to increase the apparent surface single scattering albedo and to decrease 
the surface backscattering fraction value. 
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Table 7: Retrieved Hapke's parameters and their standard deviation (the underconstrained parameters are indicated by 
"(-) ") from CRISM measurement at 750 nm (present results) compared to those from HRSC measurements at 675 nm ( Jehl 
et al.\ 20081. The yellow unit is associated to the Spirit landing site. The case 1 is the Hapke's parameters determination 
when the opposition effect parameters are set free and the case 2, when the opposition effect is neglected, with phase 



angles larger than 20° (see Jehl et at. ( 2008 1 for more details) 



Instrument 


ROI or unit 


UJ 


b 


c 


e 

(deg.) 


Bo 


h 




ROI I 


0.71 (0.05) 


0.22 (0.18) 


0.54 (0.18) 


15.78 (2.65) 


0.49 (-) 


0.54 (-) 


CRISM 


ROI II 


0.72 (0.05) 


0.27 (0.15) 


0.56 (0.16) 


15.62 (2.43) 


0.49 (-) 


0.52 (-) 






ROI III 


0.69 (0.04) 


0.19 (0.14) 


0.66 (0.18) 


13.96 (4.40) 


0.50 (-) 


0.53 (-) 




ROI IV 


0.79 (0.07) 


0.59 (0.27) 


0.56 (0.16) 


10.88 (5.11) 


0.50 (-) 


0.48 (-) 


HRSC 


case 1 
case 2 


0.72±0.02 
0.80±0.02 


0.06±0.02 
0.22±0.04 


0.34±0.06 
0.41±0.06 


18.5±1.5 
17.2±3.0 


0.73±0.07 


0.75±0.13 



4.2 Validation of Meridiani Planum 

4.2.1 Comparison to experimental measurements on artificial particles 

Our investigation shows that parameters h and c estimated from the inversion of both individual and combined (cf. Figure 
12 1 CRISM FRTs are consistent with the L-shape defined from laboratory studies. However, the bidirectional reflectance 
sampling is not adequate to provide accurate estimation of the scattering parameters, even by combining up to three 
FRTs. In this case, according to parameters h and c, the surface represented by ROI I is slightly forward scattering with 
a relatively broad lobe (low 6 and c) which is consistent with artificial materials composed of round and clear sphere or 
agglutinates and close to irregular or round rough, opaque and solid natural particles (cf. Figure 12 1. 



4.2.2 Comparison to in situ measurements taken by Pancam/MER-Opportunity 

Similar to the work done at the Spirit's landing site, a photometric study was done at Meridiani Planum by using the 
Pancam instrument from Eagle crater to Purgatory ripple in ( Johnson et al.\ 2006b). Units were defined, for example, 
as "Spherule soil" (typical soil with abundant spherules), "Bounce marks" (soil compressed by the airbags at the landing 
site), and "Outcrop rock" (bedrock). Our results are compared to two photometric sequences taken by Pancam that are 
near the selected ROI: South of Voyager (Sol 437-439) and Purgatory region (Sol 449-473) (cf. Figure |2|. In addition 
to the Spherule soil unit, Johnson et al. ( |2006b I defined the soil class "Striped soil", which characterizes soils with a 
striped appearance on the faces of some dune forms. Especially for the Purgatory region, large ripples were present and 
a specific class was defined. Bright soil deposits among the Striped soil unit were defined as a Dusty soil unit which was 
limited in spatial extent. Results are summarized in Table [8] In this article the photometric results are obtained from the 
combination of four FRT observations using the H93 version of the Hapke's model. In this way, we use the same model 



as Johnson et al. (2006b I 



Table 8: Retrieved Hapke's parameters and their standard deviation (the unconstrained parameters are indicated by 
"(+,-)") from the 2-term HG model for Sol 437-439 (Spherule Soil and Striped soil) and Sol 449-473 units (Spherule Soil, 



Site 


Unit 




S (deg.) 


b 


c 


Number 


a (c 


eg.) 


South of Voyager 


Splierule soil 


0.53 (+0.02, -0.07) 


14 (^ 


-1.-2) 


0.249 (+0.022, -0.031) 


0.491 (+0.058, -0.057) 


144 


~5- 


140 


(Sol 437-439) 


Striped soil 


0.56 (+0.01, -0.01) 


15 (^ 


-1. -1) 


0.305 (+0.010, -0.022) 


0.353 (+0.033, -0.032) 


119 


~5- 


140 


Purgatory region 


Spherule soil 


0.51 (+0.00, -0.01) 


10 (^ 


-0,-1) 


0.230 (^0.006, -0.007) 


0.761 (+0.023, -0.013) 


686 


~0- 


135 


(Sol 449-473) 


Striped soil 


0.52 (+0.24, -0.00) 


11 (^ 


-0, -1) 


0.117 (-0.011, -0.012) 


1.000 ( + , -) 


104 


~0- 


135 




Dusty soil 


0.66 (+0.02,-0.02) 


16 (- 


fl,-l) 


0.449 (+0.011,-0.020) 


0.364 (+0.019,-0.053) 


64 


~0- 


135 




Ripples 


0.52 (+0.01,-0.01) 


9 h 


-l.-l) 


0.255 (-0.008,-0.007) 


0.434 (+0.023,-0.014) 


234 


~0- 


135 



• Single scattering albedo: Figure 13 plots means and uncertainties as well as the PDF of the single scattering albedo 



UJ of the selected ROI against the values of the Spherule Soil and Striped Soil units defined by Johnson et al. { 2006b I 
at South of Voyager and the additional Dusty soil and Ripples units in the Purgatory region. The retrieved values 
of u! from CRISM are higher than those estimated from Pancam measurements except for the spatially sparse Dusty 
soil unit (first maximum of the PDF around 0.65). 



• Phase function: Figure 14 plots the asymmetric parameter b and the backscattering fraction c solutions of the 



Bayesian inversion for the ROI I in comparison to results of two photometric units (i.e.. Spherule soil and Striped 
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Figure 12: Same as Figure |8] for Meridiani Planum using observations FRT95B8, FRT334D and FRTB6B5 and their 
combination (ROI I). 
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Figure 13: Mean and uncertainties of the single scattering albedo to estimated from Pancam measurements at 753 nm for 
different geological units at South of Voyager (Sol 437-439) and Purgatory region (Sol 449-473) ^Johnson et al.\ |2006b[ ) 
compared to those estimated from CRISM measurements at 750 nm (note that the error bar represents 2a uncertainties) 
derived from 2-term HG models. The PDFs of the parameter oj estimated from the last 500 iterations of the bayesian 
inversion are also represented on the left side. This is helpful when mean and uncertainties are not representative of the 
PDF (such as the PDF of the parameter ui for the ROI I which shows a bimodal distribution). 
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Figure 14: Probability density map of the asymmetric parameter b (horizontal axis) and the backscattering fraction c 
(vertical axis) solutions from the last 500 iterations of the Bayesian inversion estimated at 750 nm obtained for the selected 
ROI. The grid is divided into 24 (vertical) x 20 (horizontal) square bins. The coloring gives the probability corresponding 
to each bin. Means and 2a uncertainties (red rhombus) are plotted too. The inversion solutions are plotted against 
the experimental b and c values pertaining to artificial particles measured by McGuire and Hapke\ ( 1995) and to natural 
particles measured by Souchon et al. ( 2011[ ). The Spherule Soil (circle), Striped Soil (square), Ripple (triangle) and Dusty 
Soil (triangle) units from Pancam sequences at South of Voyager (green and Purgatory region (blue) estimated at 753 nm 
( Johnson et al.\ 2006b I are plotted here. 



soil) from South of Voyager and three photometric units (Spherule soil. Dusty Soil and Ripples) for the Purgatory 
region. Note that the parameter c is underconstrained in the Striped Soil unit in the Purgatory region area because 
of the relative lack of phase angle coverage (~45-110°). As mentioned in Subsection |3.2[ parameters b and c axe not 
well constrained in the CRISM case but exhibit slightly more forward scattering properties than at Gusev crater. 
This results is still potentially consistent with all units in Purgatory region and South of Voyager, except for the 
Spherule soil which exhibits more backscattering properties. 

Macroscopic roughness: Means and uncertainties as well as the PDF of the macroscopic roughness 6 modeled from 



CRISM data is higher than that retrieved at the South of Voyager and at Purgatory region (cf. Figure 15 1. 



Two principal points from the previous comparison can be outlined. First, at the effective spatial scale achieved 
from space by CRISM for the photometric characterization (460 m), the surface behaves particularly like the Dusty Soil 
unit defined by Johnson et al. (2006b I at the Purgatory region: higher single scattering albedo, more forward scattering 



materials, and rougher surfaces. Second the photometric properties estimated from CRISM data are not strongly supported 
by those determined from Pancam. In fact, the selected ROI and the Pancam photometric sequences (South of Voyager and 
Purgatory region) are located several kilometers to the south of the Opportunity landing site where a geological transition 
region is observed. Indeed bright plains which exposed a smaller areal abundance of hematite, brighter fine-grained dust 
rich in nanophase iron oxides (such as in Gusev Crater plains) and a larger areal of bright outcrops are observed from 
orbital measurements (cf. Auxiliary Materials - Figure 1) compared to plains near the landing site (Arvidson et al. 2006). 



Inside the ROI I, even if the bright material unit does not seem to be the main geological unit at centimeter spatial scale. 
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Figure 15: Mean and uncertainties of the macroscopic roughness 9 estimated from Pancam measurements at 753 nm for 
different geological units at South of Voyager (Sol 437-439) and Purgatory region (Sol 449-473) ^Johnson et al.\ |2006b[ ) 
compared to those estimated from CRISM measurements at 750 nm (note that the error bar represents 2a uncertainties) 
derived from 2-term HG models. The PDFs of the parameter 6 estimated from the last 500 iterations of the Bayesian 
inversion are also represented on the left side. This is helpful when mean and uncertainties are not entirely representative 
of the PDF. 



it seems to be more abundant from the orbit. 



5 Discussion 

The previous results show that the surface photometric parameters retrieved in this study are consistent with those 
derived from in situ measurements. This cross-validation demonstrates that MARS-ReCO is able to estimate accurately 
the surface bidirectional reflectance of Mars. Nevertheless one could question the significance of these estimates. Therefore, 
this section aims at testing the non-Lambertian surface hypothesis used by MARS-ReCO and the influences of the surface 
bidirectional reflectance sampling on the determination of the photometric parameters. 



5.1 Non-Lambertian surface hypothesis 

The compensation for aerosol contribution represents a great challenge for photometric studies. In planetary remote 
sensing, the Lambertian surface assumption is often adopted by the radiative atmospheric calculations that allow retrieval 
of surface reflectance ( Vincendon et al. 2007 McGuire et al. 20081 Brown and Wolff 2009 Wiseman et al. 20121. The 



surface reflectance is then assumed to be independent on geometry (i.e., the variation of TOA radiance with geometry 
is supposed to exclusively relate to aerosol properties). This hypothesis is generally used as it simplifies the radiative 
transfer modeling. However, it has been proved that most surface materials (e.g., minerals and ices) have an anisotropic 
non-Lambertian scattering behaviors {de Grenier and Pinet 1995 Pinet and Rosemherg 2001 Johnson et al. 2006a|b 



Johnson et al. 2008 Lyapustin et al. 2010). Consequently, considering a Lambertian hypothesis can potentially create 



biases in the determination of the surface reflectance. 



For the AOT estimate some assumptions regarding the surface properties were taken in Wolff et al. (2009 1 's work 



Indeed this method assumes a non-Lambertian surface to estimate the AOT by using a set of surface photometric pa- 
rameters \ Johnson et al. 2006a|b I that appears to describe the surface phase function adequately for both MER landing 
sites (roughly associated to the brighter or dusty soils). This assumption is qualitative but reasonable for several reasons 
enumerated by ' Wolff et al. (20091. However bias in the estimation of surface bidirectional reflectance can appear. 



For 



both MER landing sites it seems from ^ Wolff et al. ( 2009| )'s work that the AOT retrievals are overall consistent with optical 



depths returned by the Pancam instrument (available via PDS). 
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Figure 16: Photometric curves corresponding to ROI I of the FRT3192 observation and composed of ten angular configura- 
tions. The reflectance values (in BRF units) are extracted at the TOA (green asterisk), and at the surface after correction 
for atmospheric contribution. The latter takes into account a non-Lambertian surface hypothesis using MARS-ReCO (red 
crosses) and a Lambertian surface assumption (blue plus sign) (see details in Subsection 5.11. In the three cases we plot 
error bars corresponding to the la uncertainty. 



Table 9: Retrieved Hapke's parameters (w, 6, c, in degree, Bq, h) and their standard deviation (the underconstrained 
parameters are indicated by "(-)") from CRISM measurement at 750 nm using our Bayesian inversion assuming (i) a 



non-Lambertian assumption (MARS-ReCO) { ^Ceamanos et al. 2013) and (ii) a Lambertian hypothesis <\Doute, ^.j [2009 1 



for the surface. The Hapke's photometric parameters corresponding to the ROI I are retrieved by Bayesian inversion as 
well as the non-uniform criterion k. {nb: number of angular configurations, g: phase angle range in degree) 



ROI 


FRT 


AOT 




6 


c 
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Bo 
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nb 


a 








A; 


k 
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k 


k 






ROI I 


3192 (non-lamb, hypo.) 


0.33±0.04 


0.68 (0.06) 
1.00 


0.17 (0.20) 
1.92 


0.62 (0.20) 
0.96 


11.62 (3.98) 
1.00 


0.52 (-) 
0.24 


0.52 (-) 
0.12 


10 


~56-112 


3192 (lamb, hypo.) 


0.33±0.04 


0.74 (0.08) 
0.99 


0.46 (0.32) 
1.12 


0.64 (0.17) 
0.74 


12.43 (4.86) 
0.98 


0.52 (-) 
0.25 


0.50 (-) 
0.03 


10 


~56-112 



In order to test whether or not the non-adoption of a Lambertian hypothesis improves the quality of the retrieved 
surface reflectance, we compare the surface reflectance and the surface photometric parameters retrieved by a radiative 



transfer-based atmospheric correction technique that adopts a Lambertian assumption for the surface ( Doute, S. 
and those given by MARS-ReCO 



Ceamanos et al. 20131 



20091 



Figure [16] presents the initial TOA reflectance and the surface reflectance derived from both atmospheric correction 
techniques as a function of the phase angle. As can be seen, the main difference is observed at high phase angle (greater 
than 90°) where lower BRF are estimated in the case of the Lambertian method. This difference can be explained by 
the non-Lambertian assumption of the surface by MARS-ReCO. Table [9] presents the photometric parameters retrieved 
from the surface bidirectional reflectance estimated from both techniques. First, solutions are found for the parameters 
to, b, c and with higher standard deviations from the Lambertian method. Second, higher b is observed from the 
Lambertian method which is not compatible with the in situ photometric results. This shows that the consideration of 
a non-Lambertian surface in the correction for aerosol contribution by MARS-ReCO is needed for the determination of 
accurate surface bidirectional reflectance and thus for the estimation of accurate surface photometric parameters. 
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5.2 Influences of the surface bidirectional reflectance sampling on the determination of 
the photometric parameters 

The accuracy of the determination of Hapke's parameters uj, b, c, 9, Bq and h highly depends on the quahty and the 
representativeness of the surface bidirectional reflectance used for their estimation. Indeed, for a given phase angle range, 
the estimation of all parameters is essentially controlled by the degree to which the angular space (incidence, emergence 
azimuth and consequently phase angles) is covered rather than by the number of available angular configurations ( Souchon 
et al. 2011 1. Also, a large phase angle domain is required. In fact, certain photometric parameter are sensitive to the 
availability of low and high phase angle values. First, the opposition effect is visible only with phase angle values less 
than 20°. Second, for an accurate determination of the macroscopic roughness, observations that span from small phase 
angle out to phase angle above 90° is needed {Helfenstein 19881. Finally, materials with strong backscattering or forward 



scattering properties require the consideration of large phase angles, ideally greater than 140° {Kamei and Nakamura 



2002, Shkuratov et al. 20071 Shepard and Helfenstein 2011 Helfenstein et al. 19911. 



As presented previously (Subsection 2.1.11, CRISM targeted observations are composed of eleven images which are 
characterized by a constant incidence angle and eleven different emergence angles. Two modes of relative azimuth are 
used corresponding to the inbound and outbound portions of the spacecraft trajectory. Based on laboratory experiments 



performed under a given incidence angle, Souchon (20121 showed that the photometric parameters are determined with a 



good reliability under the condition of varied emergence and azimuth angles, resulting in varied phase angle values. How- 
ever, each CRISM observation has intrinsically limited azimuthal and emergence coverage, which affects the determination 
of all the photometric parameters such as the opposition effect parameters, Bq and h, the phase function parameters, b and 
c, and the macroscopic roughness, 0. The solution used in the present study is to combine several targeted observations 
with different illumination conditions taken at different times in a year (because of the MRO sun-synchronous orbit) in 
order to enrich as much as possible the phase angle range by adding low and high phase angles. Nevertheless, low (less 
than 20°) and high (greater than 120°) phase angle values are not available here. 

Thanks to Bayesian inversion, the shape of the a posteriori PDF is known and will inform us whether the bidirectional 
reflectance sampling is sufficient to estimate accurate photometric parameters. In fact, a uniform or bimodal PDF are the 
consequence of a lack of the bidirectional reflectance sampling as conflrmed by synthetic tests. 

Note that the surface bidirectional sampling also influences the AOT estimate notably through the determination of the 
single scattering albedo. Large AOT error bars are related to CRISM targeted observations presenting limited geometries 
Wolff et al.\ ,2009) . However, the surface bidirectional reflectance estimated by MARS-ReCO is sensitive to the accuracy 



of the AOT input ( Ceamanos et al. 20131. 



5.3 Comparison of surface photometric results from Pancam versus CRISM 

The differences regarding the acquisition geometries between Pancam and CRISM can influence the determination of the 
photometric properties. Ideally, in order to properly compare both results, the bidirectional reflectance sampling should 
be the same. However, this cannot be achieved because of technical constraints on the acquisition of the bidirectional 
reflectance measurements in the case of a spaceborne instrument and an in situ vehicle, respectively. In addition, the 
atmospheric contribution does not affect the observation in the same manner. In the orbital case the light travels twice 
across the atmosphere while in the in situ case, it only travels once before reaching the sensor. 

Comparison of orbital photometric results to in situ photometric results demonstrates two main differences arising 
from different acquisition modes. The flrst comes from the sampling of the bidirectional reflectance. Pancam acquires a 
photometric sequence by varying the illumination angles (observations acquired at several times of day) and by varying 
the local azimuth and emission angles among similar units to order to enrich the phase angle range. 

The second difference comes from the spatial scale that is accessible. CRISM photometric curves are obtained at the 
hectometric spatial scale on an extended area whereas Pancam observes at the centimeter spatial scale. Some characteristics 
of the surface can be observed at a hectometric scale but not at centimeter scale and conversely. From the ground, Pancam 
is able to distinguish rocks and soils with relatively high variability, whereas CRISM observes extended areas for which the 
variability of rocks and soils may appear relatively reduced. Despite these differences, the choice is made here to consider 
the Pancam in situ photometric results as the ground truth, taken as a reference. 



5.4 Interpretation of the retrieved photometric parameters from orbit 

Consistent trends were found with regard to the in situ photometric results. At the hectometer-scale, the surface behaves 
like the Soil unit deflned by Johnson et al. (2006a I. The signal acquired by CRISM mixes several contributions coming 
from distinct surface materials occurring at the sub-pixel level. This mixing can be (i) linear or (ii) non-linear. In the case 
of linear mixing (geographic mixing), the measured signal is a linear sum of all the radiated energy curves of materials 
making up the pixel. Consequently, the consistency of our results with the Soil unit behavior can be due to their dominant 
relative abundance within the CRISM pixel. In the case of non-linear mixing, the signal measured is the result of a 
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non-linearly weighted combination of signals coming from different units. In both cases, the apparent consistency with 
the Soil unit could be coincidental. 

Previous orbital observations acquired at the MER landing sites gave several clues concerning its physical and chemical 
properties. This information can be used to validate the present photometric results to provide a better understanding of 
the geological meaning of the estimated photometric properties. First, we compare our results to the albedo estimates by 
Thermal Emission Spectrometer (TES) instrument. Spirit landed in a relatively dark streak feature which is composed 
of hundreds of dark-toned, small, sub-parallel streaks created by dust devils, roughly oriented NW-SE. This feature is 



et al 2004 



characterized by low albedo values (between 0.16 and 0.22) ( Martinez- Alonso et al. 20051. Pancam in situ analysis (Bell 
Farrand et al. 2006[ ) and coupled OMEGA measurements and Spirit data ( Lichtenberg et al. 2007 ) show 



that the Spirit landing site is mainly composed of soil deposits sprinkled with small rocks interspersed with dust deposits 
composed of fine-grained nanophase iron-oxide-rich. These results are consistent with the intermediate single scattering 
albedo estimates modeled from CRISM observations {uj ^ 0.68). Meridiani Planum is characterized by a low albedo 
estimated at 0.09-0.19 [Mellon et al. 2000). Dust-free dark sand surface is mainly observed on the Meridiani Planum 



landing site. However the intermediate single scattering albedo value modeled from CRISM dataset is more consistent 
with bright plains which exposed a smaller areal abundance of hematite, brighter fine-grained dust rich in nanophase iron 
oxides (such as in Gusev Crater plains) and a larger areal of bright outcrops are observed several kilometers to the south 



of the Opportunity landing site from orbital measurements ( Arvidson et al. 2006 1 



Second, we compare our results to the thermal inertia (TI) estimated by THermal EMission Imaging System (THEMIS) 
images. This is a measure of the resistivity of surface materials to a change in temperature which is function of the particle 
size, bulk density and cohesion. The intermediate THEMIS-derived TI at Gusev landing site (186-347 SI, 3x6 km resolution, 
[Martinez- Alonso et al. 2005| )) suggests that the surface is dominated by unconsolidated materials (cohesionless surface) 
which is consistent with the photometric estimates. This shows that the surface dominantly behaves like the Soil unit, 
particularly those observed at the Landing site which exhibit more small clasts in soil. At Meridiani landing site, lower 
THEMIS-derived TI are found (i.e., 175 SI, 3x6 km resolution, ( Jakosky et al. 2006 1) suggests that the surface is dominated 
by more fine-grained materials. Those pieces of information are consistent with the CRISM-derived photometric estimates 
showing that the area has similar properties as the Soil unit. 

Even if consistent trends were found, one difference is observed with regard to the in situ photometric results. The 
difference concerns the macroscopic roughness 0. At Gusev Crater, the intermediate values of conjugated to the presence 



of more small clasts on soil ( Ward et al. 



2005 



Johnson et al. 2006a I suggests that the parameter 9 is mainly sensitive to 



micro-structure. In contrast, slightly rougher soil was modeled from CRISM data at Meridiani Planum. This outcome can 
be explained by (i) the presence of rocks but very few rocks are observed ( Golombek et al. 2005 1 and/or (ii) the presence 



of less than 1 cm-hematite-rich concretions ( Calvin et al. 2008 1 and/or (iii) the presence of morphological structures such 
as cm- high ripples ^Golombek et al. 2010). The second and third explanations could be responsible for the roughness 
measured at hectometer scale, suggesting that the parameter 9 is sensitive to both the specific local micro- and macro- 
structure. At the Gusev plain, it seems that the surface roughness is dominated by small clasts in soils, whereas, at 
Meridiani Planum it is dominated by the soil particles (concretions) and ripples. 



6 Conclusions 

CRISM observations acquired over the MER-Spirit and MER-Opportunity landing sites allow us to validate the accuracy of 



(i) the Martian surface bidirectional refiectance estimated by the procedure named MARS-ReCO developed by Ceamanos 



et al. (2013 1 and (ii) the determined Martian photometric parameters, by comparing our results to the in situ photometric 
results modeled from data taken by the Pancam instrument on-board both rovers. Indeed, Hapke's parameters (a;, 9, b and 
c) estimated from CRISM measurements (using a combination of CRISM multi-angle images after individually by applying 
to each FRT observation the MARS-ReCO atmospheric correction procedure), are mainly consistent with parameters 
modeled from in situ measurements taken at the Spirit and Opportunity landing sites. The innovative assumption of a 
non-Lambertian surface is used in our methodology to accurately estimate intrinsic surface photometric properties from 
space. Our results appear to improve (i) those achieved in HRSC-based photometric studies, in which no aerosol correction 
was used, and (ii) those estimated from surface bidirectional refiectance derived from atmospheric correction assuming a 
Lambertian surface, as they better compare with the in situ Pancam results. This outcome of our study shows that MARS- 
ReCO gives access to consistent surface bidirectional refiectance. As a consequence, surface photometric parameters can 
be reliably estimated from CRISM observations, provided that the atmospheric conditions are not turbid. As presented 



previously, the AOT values derived from Wolff et al. (2009)'s work are estimated with assumptions concerning the surface 
properties, imposed to be similar as MER observations at both landing sites. Further developments could lead to the joint 
estimation of AOT and surface bidirectional refiectance by considering the full spectral dimension. 

This presented approach may suffer from intrinsic limitations due to the scarcity of CRISM measurements with broad 
bidirectional refiectance sampling. However the combination of several CRISM observations alleviates this problem by 
improving the phase angle range thus better constraining the determination of photometric parameters. Furthermore, 
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Since September, 2010, the inbound segment in targeted observations is absent due to problems of the instrument gimbal 
( Murchie 2012 1. Only 6 angular configurations are still available after this date. Consequently, more targeted observations 
must be combined in order to compensate this limitation. 

The presented methodology opens the possibility to map the surface bidirectional reflectance thus the spatial variations 
of the photometric parameters. The determination of the physical state of the surface materials (i.e., mean grain size, grain 
sphericity and grain heterogeneity, mean surface porosity and surface brightness) through the study of their photometric 
properties (phase function of the particle, single scattering albedo, surface roughness, and opposition effects) provides 
complementary information to composition. It has implications for the characterization of different types of terrain over 
Mars (i.e., extrusive and intrusive volcanic terrains, sedimentary terrains, impact craters) and for identification of surface 
temporal modifications (i.e., aeolian process, space weathering process). Moreover, the present work has implications 
for spectroscopic interpretation as it gives information about the particle (i.e., size) and grain organization which have 
significant infiuence on spectral properties (Poulet and Erard 2004). 

In future work, both challenges of map creation and improvements for non-flat terrains shall be taken into account, 
with a particular emphasis on the complex handling of local slope effects which is not straightforward in the case of 
multiple reflection between surface facets. Note that MARS-ReCO may be used for mapping the photometric properties 
of the Martian surface even when using a single CRISM targeted observation. In this case, however, the photometric 
inversion should be focused to retrieve only single scattering albedo for the best surface bidirectional reflectance sampling. 
Furthermore, use of a single-lobed Henyey-Greenstein function would reduce the number of unknown photometric param- 
eters and therefore could allows us to constrain the asymmetry parameter besides the single scattering albedo. Likewise, 
laboratory studies bearing on the photometric behavior of natural surfaces under controlled conditions and under a large 



range of observational geometries {Souchon et al. 20111 are necessary to better document and interpret the variations of 



the photometric parameters in relation with surface geological processes such as eolian, aqueous and impact processes. 



Acknowledgments 

This work was supported by the French Space Agency CNES (Centre National d'Etudes Spatiales) and PNP (Programme 
National de Planetologie) from INSU (Institut National des Sciences de I'Univers). The authors would like to thank 
Michael Wolff for making his aerosol optical thickness values available for this study. We would like to gratefully thank 
Jeffrey Johnson, Mathieu Vincendon and the anonymous reviewer for their constructive comments that substantially 
improved this article. 



Appendix A 

Non-uniform PDF Central moments /i„ (such as the variance /i2 of order two) are commonly used for statistical 
purpose while cumulants fc„ have the advantage to present unbiased statistical estimator for all orders ( Fisher\ 1930). 



Also, the flrst three cumulants are equivalent to the central moments. For a random variable following the PDF f{x) in 
(0,1), the cumulant generating function is: 

1 

$(t) = lnTi^{/(a;)} = In j e"^/(a:)dx (9) 



Cumulants of order n are deflned by : 



kn = ^+o[t) (10) 



The first four cumulants of a uniform PDF are: 

Ml (11) 

= (12) 

M3 (13) 

fc4 = -Y^=M4-3(/i2)' (14) 
Thus we propose to estimate the non-uniformity of the results of the inversion with: 
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fci-12 fc2-112 ^3 fc4 + 1120 
' 12 ' 112 '160 1120 ' ^ ' 

We perform 10,000 uniform random vectors of 500 samples (identical to the inversion procedure). Since the maximum 
is fc=0.47 for the most extreme event, we propose to have non-uniform PDF for A;>0.5. For the inversion purpose, since 
the a priori PDF's on the parameters are uniform if the results of the inversion on one parameter has fc<0.5, we conclude 
that this parameters is not constrained by the observations. 
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